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Traveling waves are ubiquitous in nature, and control the speed of 
many important dynamical processes, including chemical reactions, 
epidemic outbreaks and biological evolution. Despite their funda- 
mental role in complex systems, traveling waves remain elusive be- 
cause they are often dominated by rare fluctuations in the wave 
tip, which have defied any rigorous analysis so far. Here, we show 
that by adjusting nonlinear model details, noisy traveling waves can 
be solved exactly. The moment equations of these tuned models 
are closed and have a simple analytical structure resembling the de- 
terministic approximation supplemented by a non-local cutoff term. 
The peculiar form of the cutoff shapes the noisy edge of traveling 
waves and is critical for the correct prediction of the wave speed and 
its fluctuations. Our approach is illustrated and benchmarked using 
the example of fitness waves arising in simple models of microbial 
evolution, which are highly sensitive to number fluctuations. We 
demonstrate explicitly how these models can be tuned to account 
for finite population sizes, and determine how quickly populations 
adapt as a function of population size and mutation rates. More 
generally, our method is shown to apply to a broad class of models, 
in which number fluctuations are generated by branching processes. 
Due to this versatility, the method of model tuning may serve as a 
promising route towards unraveling universal properties of complex 
discrete particle systems. 

The wave-like spread of discrete entities pervades our ev- 
eryday life. For example, the spread of ions, pathogens 
and beneficial mutations control the human heart beat, the 
yearly threat of influenza and evolutionary progress [T]. A 
thorough understanding of how these waves form and spread 
has numerous applications ranging from the control of chemi- 
cal reactions 2 , 3 to the prediction of epidemic outbreaks |4] 
[5] . Great research effort has therefore been made on the ques- 
tion of how traveling waves emerge in complex systems from 
the multiplicity of relatively simple interactions, in particular 
the random dispersal of particles and reactions between par- 
ticles. To simplify the analysis, most theoretical studies have 
been neglecting number fluctuations, which are inevitable in 
systems of discrete particles. However, when stochastic sim- 
ulations became feasible, it was found that those previously 
ignored fluctuations can have a strong impact on the dynam- 
ics of waves [H [T] [8] . Simple models of biological evolution 
(described below) serve as the prime example for this dras- 
tic sensitivity on noise because they are dominated by the 
few most fit individuals in a population, and break down if 
number fluctuations are neglected jGj. With the added diffi- 
culty of particle discreteness, the analysis of traveling waves 
became one of the important challenges of statistical physics 
and mathematical biology, which still defies systematic ana- 
lytical techniques. 

Here, we show that an exact analysis of noisy traveling 
waves is feasible when the nonlinear details of the underlying 
model are chosen appropriately. The main idea is to tune the 
nonlinearities to obtain the least difficult math, while retain- 
ing the universal features of the model. Specifically, we show 
that it is possible to close the hierarchy of moment equations 
using a suitably designed nonlinear constraint on the dynam- 
ics. The tuned model can then be used to extract the universal 
features of noisy traveling waves, for which only heuristic ap- 
proaches have been available so far. The method of model 
tuning can be naturally generalized to a wide range of un- 
solved stochastic nonlinear problems, and therefore provides 



a promising new tool to unraveling universal features of non- 
equilibrium systems. 



Noisy wave models 

Simple models of evolution [6l [H [TO] [TTJ [12] [13] provide spec- 
tacular examples of noisy traveling waves because of their 
drastic sensitivity to rare fluctuations in the wave tip. As 
illustrated in figure [l] these waves describe the continual in- 
crease of growth rate, also called fitness, due to spontaneous 
mutations in a finite population. The wave speed, which is a 
measure for how quickly populations adapt, increases without 
bound if number fluctuations are neglected [6]. To reproduce 
a steady state with constant wave speed observed in stochastic 
simulations, any analysis has to incorporate number fluctua- 
tions, at least heuristically by introducing an ad hoc "cutoff" 
in the noisy tip of the wave 6 . The relationship between 
traveling speed and population size has been investigated ex- 
tensively in the recent literature 14, 15 , 13, 16 , with some 
controversy regarding the universal behavior 17 , IS^ , in order 
to interpret growth rate measurements in microbial evolution 
experiments [T9l[20] . 

We now take a closer look at the models underlying these 
fitness waves to elucidate their characteristic structure. This 




Growth rates (fitness) 

Fig. 1. A paradigmatic example for noisy traveling waves are "fitness" waves 
arising in simple models of evolution. The colored particles represent individuals with 
characteristic growth rates, or fitnesses (horizontal axis). Individuals can mutate, 
replicate ("birth") and be eliminated from the gene pool ("death"), as illustrated. 
These simple dynamical rules give rise to a distribution of growth rates resembling a 
bell-like curve at steady state, which propagates towards higher growth rates like a 
solitary wave. The random fluctuations in the tails of the wave have precluded any 
rigorous analysis in the past. 



Reserved for Publication Footnotes 



www.pnas.org 



PNAS I Issue Date | Volume | Issue Number | 1-6 



will lead us to a general model of noisy traveling waves, which 
forms the basis of our analytic approach. Most evolution mod- 
els implement the reproduction of individuals by a standard 
branching process. The growth rate is subject to small vari- 
ations due to random mutations, which are modeled through 
a random walk, for instance by standard diffusion [6]. In ad- 
dition, all evolution models contain a nonlinearity that limits 
the total population size. This accounts for the fact that nat- 
ural populations must stay finite due to resource limitations. 

The generic mathematical structure of such evolution 
models can be framed as follows. The state of the popula- 
tion at time t is described by a function Ct{x) that represents 
the number density of individuals having a growth rate x. The 
population is assumed to evolve in discrete timesteps of size 
e, which is eventually sent to zero in order to obtain a contin- 
uous time Markov process. Any such timestep consists of two 
sub-steps. The first one, illustrated in Fig. [2}\, embodies the 
mutational process, which leads to slight variations in growth 
rates, and the process of reproduction, which causes individ- 
uals with growth rate above the mean to increase in relative 
number. In general, the combined effect of both processes can 
be described by the equation 

ct+e - ct = eCct + V^ectT] , [1] 

which takes the concentration field from ct to an intermediate 
value ct+e- The term eCct represents the expected change in 
density due to reproduction and mutations. This term is lin- 
ear in the concentration field because the number of offspring 
and mutants per timestep is proportional to the current pop- 
ulation density. The linear operator C is similar to a Hamilto- 
nian in physics. A specific example will be discussed later on. 



1. Sub-Step: Reproduction & Mutations 




Growth rates X 



Fig. 2. A computational timestep in fitness wave models consists of two sub-steps. 
The first sub-step (upper panel) is reproduction and mutations. Growth favors the 
most fit individuals as indicated by the red arrows. The resulting population density 
c(x) will generally violate the constrained of constant population size. The second 
sub-step (lower panel) restores the prescribed population size by a random elimination 
of individuals from the population (this step may also be interpreted as Darwinian se- 
lection). The result of the whole timestep is a constrained branching random walk of 
"genotypes", which shifts the fitness distribution towards higher fitness. 



The stochastic term \/2ectr] in equation accounts for all 
random factors that influence the reproduction process. The 
function r]{x) represents standard white noise, i.e., a field of 
delta correlated random numbers, r](x)r](y) — 6(x — y), where 
/ denotes the ensemble average of a random variable /. The 
amplitude oc ^/ect of the noise term in equation ([T]) is typical 
for number fluctuations: Due to the law of large number, the 
expected variance from one timestep to the next is propor- 
tional to the number ect of expected births or deaths during 
one timestep. 

Because the reproduction step changes particle numbers, 
another sub-step of population control is required to enforce 
a constant population size. In most models and experi- 
ments [191 l2Qj . this step is realized by a random culling of 
the population: individuals are eliminated at random from 
the population until the population size constraint is restored, 
see Fig. [2^3. Mathematically, the selection step can be cast 
into the form 

Ct+e=Ct+e(l-A) , [2] 

where A represents the fraction of the population that has to 
be removed to comply with the population size constraint. 
The second sub-step completes the computational timestep, 
and takes the concentration field from the intermediate state 
Ct+e to the properly constrained state Q+e- 

The combination of noise and nonlinearity (via the popu- 
lation constraint) has led many theoretical studies to engage 
in approximations that are often hard to justify or to con- 
trojj Here, we take a different approach by optimizing the 
form of the nonlinearities in order to obtain an exactly solv- 
able model. To this end, we first generalize the above model 
in the following way. Instead of enforcing a fixed population 
size, we allow for a whole class of constraints. Specifically, we 
assume that the selection step enforces a constant value of 1 
for the inner product of the concentration field Ct(x) and a 
new function u{x), 



Here, we have introduced the "bra-ket" notation commonly 
used in quantum mechanics. The function u{x) defines the 
selection rule and will be called selection function from now 
on. It will serve as a "tuning wheel" in the following. If one 
chooses u = , one obviously recovers the fixed popula- 
tion size constraint. For any other choice, the population size 
will not be fixed. At best, one obtains a steady state with a 
population size fluctuating around its mean value, N. 



Results 

Model tuning. Equations ([l] |2] [3| define a class of models 
that generate branching random walks subject to a global con- 
straint. As we argue in the Discussion section, these features 
not only characterize fitness waves but furnish the essence 
of most noisy traveling waves. The fundamental difficulty in 
analyzing such models becomes apparent when we try to de- 
termine the typical dynamics by averaging over the stochastic 
noise term. This can be done by eliminating A using the con- 
straint, equation ([3|, sending e to zero and carrying out the 
ensemble average. This straight forward calculation, detailed 
in the Supporting Information, yields 



dtct = - 2u) Ct - {ct I (£t - 2u) u)ct [4] 



remarkable exception is the exactly solvable "exponential" model by Brunet et al. |21| 
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for the mean concentration field ct. Equation Q reflects the 
usual "horror" inherent to nonlinear stochastic problems: Mo- 
ment equations do not close in general. The mean depends 
on the second moment, the second on the third, and so on. A 
whole hierarchy of moment equations would have to be solved 
self-consistent ly to make progress. 

However, if we consider u{x) as a tunable function, equa- 
tion Q allows for different kind of simplification. Suppose, 
the solution u^{x) of 

(^C"^ - 2u.^ = [5] 

exists, and we choose (x) as the selection function. For this 
particular model, the nonlinearity in equation Q disappears 
identically. Thus, the dynamics of the first moment becomes 
linear, 

dtct = {C-2u,)ct , [6] 

and hence solvable. In particular, if a steady state exists, 
the steady state concentration field c is in the null space of 
C — 2u^ . These observations suggest an algorithm to construct 
for a given linear operator £, a solvable constrained branch- 
ing walk model: First, identify the selection function u^(x) 
for which the averaged dynamics becomes simple (i.e. linear). 
To this end, one has to solve equation ([5]), which is in general 
a deterministic partial differential equation. Second, solve the 
corresponding moment equation ([6|, which is guaranteed to 
be linear for the chosen selection function (x) . 

At this point it is not clear, whether equations ([5] [6| have 
stationary solutions, and, if so, whether these solutions can be 
used to extract universal features of traveling waves. There- 
fore, we apply our recipe to a simple model of fitness waves. 
This model is defined by the operator 

£evo = (X - Xo) + Vda: + DO^ , [7] 

which corresponds to a branching random walk on a reaction 
rate gradient 22 . It has been proposed as the simplest model 
of asexual evolution in order to explain the fitness growth ob- 
served in directed evolution experiments with large popula- 
tions of RNA viruses [6]. The term (x — xo)ct in equation ([7| 
represents the reproduction with xq being the mean growth 
rate of the population. The mutational process is modeled as 
a diffusion process with diffusivity D. This diffusion approxi- 
mation is justified when mutations occur at a higher rate than 
the typical growth rate difference they confer, which is an ap- 
propriate assumption for viruses with large mutation rates. 
For microbes with lower mutation rates, it is necessary to dis- 
cretize the fitness landscape and model mutations as fitness 
jumps of finite size. The discretized scenario is slightly more 
complex as it has one additional parameter, namely the char- 
acteristic fitness effect of mutations, but can be discussed in 
complete analogy to the present case of a continuous fitness 
landscape (Supporting Information). Finally, the term vdxCt 
appears in equation ([7| because the dynamics of the fitness 
wave is described in a reference frame moving with the average 
speed V of the fitness wave. 

Using £evo in equation (jsj, we find that the equation 

Ddlu^ - vdxu^ +xu^-2ul = [8] 

specifies the selection function u^{x). Once i^* is determined, 
the mean concentration field follows as the stationary solution 
of the linear equation (|6|. This second task requires no spe- 
cial effort if the random walks are simple diffusion processes 
as in equation ([7|. Then, c at steady state and u^^ are simply 
related by 

c (X. u^e ' ^ [9 J 



which generates equation (|8| when inserted into equation (|6| . 
Thus, c follows from directly by multiplying with a decay- 
ing exponential. Theprefactor in equation ^ is set by the 
constraint equation (Is]). 

The set of equations (|8]|9| is controlled by just one param- 
eter vD~^l^ ^ as can be seen by introducing scaled variables. 
In Fig. [3] we depict representative numerical solutions for 
both regimes of large and small values of this control param- 
eter, and compare with simulations of the tuned model. We 
find that both regimes exhibit strongly different wave profiles. 
Whereas the mean concentration c exhibits for the most part 
a Gaussian shape in regime of large wave speeds (vD~'^^^ > 1 
), it is strongly skewed for small speeds. Notice the nearly per- 
fect agreement between numerics and stochastic simulations, 
which confirms our theoretical framework and demonstrates 
the feasibility of our model tuning approach. 

Interpretation of the tuned models. Among the three graphs 
depicted in Fig. Is] only the function c{x) has an obvious 
interpretation as the fitness wave profile. Using the results 
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Fig. 3. Stochastic simulations confirm exact theoretical predictions. The mean 
population density c and the functions li* (x), g = u^c for two different choices for 
the parameters of the evolution model are depicted. The selection function u^{x) 
has the intuitive interpretation of a fixation probability, see main text. The distri- 
bution g = u^c indicates from which co-ordinate the individual is sampled whose 
descendants will eventually take over the front. To obtain the depicted curves, we 
first identified u^{x) by solving equation JsJ numerically and then determined c{x) 
by stochastic simulations of the tuned model with selection function w(a::) = u^{x). 
Notice the near perfect agreement between stochastic simulations (symbols) and the- 
ory (solid lines). The control parameter for A) and B) was set to vD~^l^ = 5.17 
{D = 10-11) and vD-'^/^ = .464 {D = 10'^), respectively. 



^ Note that, as required for a probability density, the integral of g(x) = u^c over x is equal to 
1 by virtue of the constraint equation jsj. 
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obtained in Ref. [23], one can also give an intuitive interpre- 
tation of the functions u^{x) and g{x) = u^c. Both functions 
relate to the phenomenon of fixation. Imagine sampling an 
individual at position x and labeling it with an inheritable la- 
bel (neutral mutation). As the dynamics proceeds, the abun- 
dance of this label will change due to number fluctuations and 
the fitness of its carriers. Eventually, this label will either go 
extinct, or become fixed in the population. The latter case 
occurs if the descendants of the labeled individual take over 
the population. 

Fixation events are much more likely if the initially la- 
beled individual belongs to the fitter part of the population. 
We thus expect the probability of fixation to be a steeply 
increasing function of x, similar to the function In- 
deed, the interpretation of u^{x) is precisely that of a fixation 
probability of a particle at position x, which is derived in the 
Supporting Information. It is interesting to note about equa- 
tion ^ determining u^{x), that a very similar equation is 
well-known to describe the survival probability of an uncon- 
strained branching random walk 12 . The only difference is 
the pre-factor of 2 instead of 1 in the nonlinearity of equation 

The product g{x) = u^c also has an intuitive interpre- 
tation in terms of a probability densitjj^ It represents the 
positional distribution function of the individual whose de- 
scendants eventually will take over the population 23 . Even 
though there must surely be such a "lucky" individual at any 
time, it cannot be described deterministically, simply because 
the fixation event depends on random future events. Thus, 
the position of the "common ancestor of future generations" 
can only be described probabilistically, similar to the position 
of quantum mechanical particle. 



Wave speed. Next, we compare the evolution model with its 
tuned counterpart, considering at first the relation between 
speed and population size. From Fig. |4] it can be seen that 
deviations are significant only for intermediate values of the 
control parameter vD~^^^ , and disappear in the limit of large 
and small values. The scaling v ^ In^^^ N for large speeds 
has been observed previously, but could only be modeled by 
introducing a cutoff into the mean- field equations ^i22 . Our 
tuned model, which exhibits the same scaling, can be used to 
rationalize and refine the heuristic cutoff idea. To this end, 
consider the closed equation for the mean population density 



at steady state. 



= DdxC + vdxC + (x — xo)c - 




In (iVDV3^ 



Fig. 4. The relation between the scaled speed, vD~^^^ , and the scaled popu- 
lation size N D^^^ as obtained numerically for the tuned model (black line), which 
is solvable, and from stochastic simulations of the original evolution model 6 with 
a rigid population size constraint (red points). Both curves approach each other in 
the limit of large and small N . The data are consistent with the asymptotic scal- 
ing V ~ In"*"/"^ N (see inset) and v N ^or large and small population sizes, 
respectively. 
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which is obtained by eliminating in equation ([8| using equa- 
tion |9| and the constraint (|3|. The closed moment equation 
( |1Q| ) has the form of the deterministic approximation except 
for the last term, which acts (for large population densities) 
similar to a cutoff: It is negligible for small x, and completely 
dominates for sufficiently large x because of the exponential 
"amplification" factor in the numerator. This leads to an 
asymptotic wave profile ^ e~^^ in the tip of the wave, just 
as found from an heuristic cutoff in the reaction rate [7]. The 
location of the cutoff can be determined by balancing the last 
term in equation (10) with the reaction term (x — Xo)c. In 
Fig. [3] this crossover point can be identified as the x- value 
where u crosses over from exponential to linear. Notice that 
the functional form of the cutoff is quite different from a step- 
function, which has been used as a cutoff in Ref. 6 and most 
subsequent studies. The peculiar nonlinear form of the cutoff 
turns out to be essential for determining the leading order cor- 
rections to the wave speed (and its fluctuations) in the limit 
of large speeds, as we show in the Supporting Information. 
In the opposite limit of small speeds ^ 0, we find that the 
tuned model exhibits v ^ N/2, which can be derived from a 
perturbation analysis of the moment equation 24 . Surpris- 
ingly, the exact same asymptotic is measured for the original 
evolution model, suggesting universal behavior not only for 
large but also for small speeds. 

Fluctuations. Another telling comparison concerns the fluctu- 
ations of waves, rather then the mean density field. At first 
sight, fluctuations in both models seem to be very different be- 
cause they arise from different statistical ensembles: Whereas 
wave speeds fluctuate in the evolution model by Tsimring et 
al. 6 , they are perfectly constant in the tuned model. On 
the other hand, population sizes are constant in the evolution 
model but fluctuate in the tuned model. One might expect, 
however, that both ensembles become equivalent in the uni- 
versal large population size limit, similar to different ensem- 
bles in the thermodynamic limit of statistical mechanics. On 
the basis of this equivalence hypothesis, one can try to infer 
the fluctuations in wave speed of the evolution model (and 
thus wave diffusivities) based on the fluctuations of the mor- 
tality rate A in the solvable tuned model (Supporting Informa- 
tion). This reproduces correctly the known diffusivity scaling 
i^wave ^ ln~^ N of Fisher-Kolmogorov waves 25 and predicts 
a novel scaling Dwave ln~^ N for evolutionary waves. We 
would like to stress, however, that these analytic results rest 
on the equivalence of the two descriptions in the large popu- 
lation size limit, which remains to be proven. 



Discussion 

As a case study for noisy traveling waves in general, our analy- 
sis concentrated on simple, yet unsolved, models of evolution. 
These models effectively describe branching random walks 
subject to a global constraint T7: The branching random 
walk has the effect of modeling the growth of the population 
and the mutations through diffusion along the fitness axis. 
The global constraint fixes the population size accounting for 
the fact that natural population must stay finite due to re- 
source limitations. 

Even these simplest evolution models could not be solved 
previously because of the nonlinearity associated with the 
rigid constraint. However, while a limit on the population 
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size is natural, there is no obvious biological reason for de- 
manding a non-fluctuating population size. After abandoning 
the idealized fixed population size constraint, we could show 
that a solvable evolution model with a closed moment equa- 
tion can be obtained when the form of the constrained is tai- 
lored to satisfy equation ([5|. The solvable model was found 
to reproduce the properties of the fixed population size model 
very well, in particular for large and small population sizes 
(or wave speeds). 

Our approach of model tuning not only applies to evo- 
lutionary waves but to any solitary wave that is generated 
by a branching random walk under a global constraint, or 
some other nonlinearity, to keep particle numbers finite. Such 
models, which we may term constrained branching random 
walks, share universal features that only depend on the form 
of the linear part of the model, as has been found in extensive 
computer simulations [8 . The universality class is defined 
by the linear part of the model, which describes a branch- 
ing process and a random walk. For instance, the important 
class of Fisher-Kolmogorov waves [261 127j has a reaction term 
that saturates at a constant value in the tip of the wave, in 
contrast to the linear reaction rate gradient that character- 
izes the above evolutionary waves. This difference results in 
a different scaling of the wave speed and its dependence on 
number fluctuations, which is correctly reproduced by our ap- 
proach beyond the leading order "cutoff" approximation 7 
(Supporting Information). Wave speeds are also strongly de- 
pendent on the stochastic particle motion, which is diffusive 
only in the simplest cases. The spread of epidemic waves in 
the modern human population, for instance, is dominated by 
the scale free air transportation network ,5_. This leads to 
super-diffusion [28], which can be incorporated in our frame- 
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work by the use of a fractional diffusion operator in £. Apart 
from these obvious extensions of our approach, our method 
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ios C{t) (Supporting Information). This allows to account for 
dynamic driving forces and to study, for instance, waves in 
temporally changing environments. 

In summary, we have developed and applied an exact 
method to determine the universal dynamics of noisy traveling 
waves. By tuning nonlinear details of traveling wave models 
without changing their universal features, we could close the 
hierarchy of moment equations for a large class of noisy trav- 
eling waves. These solvable models are amenable to a simple 
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bles mean field theory supplemented with a cutoff, and thus 
provides, to our knowledge, the first rationale for the heuristic 
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for the correct prediction of the wave speed and its fluctua- 
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finite. This encompasses simple models of microbial evolution 
and range expansions (Fisher-Kolmogorov waves), for which 
we explicitly demonstrated the utility of our approach. It will 
be interesting to see whether model tuning can be extended 
to situations where interactions are inherently local, such as 
in the problem of directed percolation 29 . 
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